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The mapping of the Nonlinear Schrodinger Equation with a random potential on the Fokker- 
Planck equation is used to calculate the localization length of its stationary states. The asymptotic 
growth rates of the moments of the wave function and its derivative for the linear Schrodinger 
Equation in a random potential are computed analytically and resummation is used to obtain the 
corresponding growth rate for the nonlinear Schrodinger equation and the localization length of the 
stationary states. 
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I. INTRODUCTION 



In this work we consider the problem of the one-dimensional Anderson localization Jl|, [^ for the nonlinear 
Schrodinger equation (NLSE) [d, U, l5|- In spite of the extensive research, many fundamental problems are still 

Q open, and, in particular, it is not clear whether in one dimension (ID) Anderson localization can survive the effects of 

D , nonlinearities. This problem is relevant for experiments in nonlinear optics, for example disordered photonic lattices 

H ■ [aj 01 5 where Anderson localization was found in presence of nonlinear effects as well as experiments on Bose-Einstein 

I \ Condensates (BEC) in disordered optical lattices | § . [3 , lOL [ill . Il3 |. The interplay between disorder and nonlinear 

C^ ■ effects leads to new interesting physics [9^10j_ 13,llJ, [lalla]- In particular, the problem of spreading of wave packets 

"^ ' and transmission are not simply related 13, [iMl, in contrast with the linear case. 

I t , We consider one-dimensional localization of stationary solutions of the nonlinear Schrodinger equation (NLSE) in a 

d ' random ^-correlated potential V{x) with a Gaussian distribution (white noise), of zero mean and variance 2D, namely, 

H , (I^(2;)I^(a;')) = 2DS(x — x'). The linear version of this model was studied extensively in the past [l^. Following our 

I ' previous analysis [2 0( 1 we study Anderson localization of stationary solutions with energies lo in the framework of the 

'^ ; stationary NLSE 

O : a;0(a:) - ^dl4>{x) + P(l>^{x) + F(x)</)(x) , (1) 

\ where (t){x) is real, the variables are chosen in dimensionless units and the Planck constant is h = 1. For the lattice 

T— I ■ version of the model, it was established rigorously [3, [2l|, [23 that the stationary solutions of this equation are 

^ ' exponentially localized for a wide range of conditions. It was also argued that the rate of growth of moments for the 

|n , stationary NLSE ([1]) coincides exactly with the linear case [20| and determines the localization length. 

^^ ' We will specifically calculate {(jP'{x)) of solutions of Eq. ([T|) that are found for a certain lo, with given boundary 

conditions at some point, for example ^(x — 0) and (j)'{x = 0), where prime means the derivative with respect to x. 

This will be done with the help of the analogy with the Langevin equation [19, 20, 23]. In particular, we will calculate 

the growth rate of the second moment. 
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§: 2,= i„t<:«i»>o, f = i, (2) 
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. , ■ that will turn out to be independent of j3, where ^ is the localization length. Note, it is different from the usually 

''^ ,. studied (in the linear case) self averaging quantity 7^ = i^(ln(/)^(x)) = ^ limjj^oo ^'''J-^'' , and 7 is a smooth 
d function of energy. Since the distribution of the random potentials is translationally invariant, it is independent of the 
choice of the initial point as x = 0. Like in the linear case, starting from a specific initial condition, (j){x) will typically 
grow. For specific values of w at some point this function will start to decay, so that a normalized eigenfunction is 
found. This is the approach a la Borland [2J, [2g| that was made rigorous for the linear case in [2fl[27|- Here, following 
[20| we extend this approach, in a heuristic form, to the nonlinear case. The envelope of the wave function will grow 
exponentially if we start either from the right or from the left. The value of uj results from the matching condition, so 
that an eigenfunction has some maximum and decays in both directions as required by the normalization condition. 
The exponential decay is an asymptotic property, while the matching is determined by the potential in the vicinity of 



the maximum. This observation is crucial for the vahdity of this approach and enables us to determine the exponential 
decay rate of states from the solution of the initial value problem ((T|) . In [20| a linear equation for the moments of 
(j){x) of Eq. ([T]) and its Langevin analog was derived and it was shown that the exponents that control the growth of 
the moments are identical to the ones of the linear system (/? = 0) . In the present paper it will be shown that this is 
also correct for the asymptotic behavior of the moments. 

The outline of the paper is as follows. The analogy with the Langevin equation and dynamics of the moments is 
outlined in Section II. The generalized Lyapunov exponents, that are the eigenvalues which determine the growth of 
the various moments, are presented in Section II and their asymptotic behavior is derived analytically for the first 
time in Appendix A (they were found numerically in [2^|). The resulting asymptotic expansion for the growth of the 
moments is presented in Section III. The results are summarized in Section IV. 

II. THE FOKKER-PLANCK EQUATION AND LYAPUNOV EXPONENTS 

Following Ref. [20|, we perform the calculation of (0^(a;)) by using the analogy with the classical Langevin equation 
[l^[23]- Therefore, considering the a:-coordinate as the formal time variable on the half axis x = t £ [0, cxd), Eq. ^ 
reduces to the Langevin equation 

+ w0 - /303 - y(r)0 = (3) 

with the S correlated Gaussian noise V{t). Now we introduce new variables u — cj) and z; = (/) = ^ (that play the role 
of position and velocity in the Langevin equation) and a distribution function of these new variables is P = P(w, v, r). 
The dynamical process in the presence of the Gaussian i5-correlated noise is described by the distribution function 
that satisfies the Fokker-Planck equation: (FPE) [13, HO] 

drP - [uju ~ f3u^]d^„P + vduP - 2Du^d^P = , (4) 

where 2Dv? is the only nonzero component of the diffusion tensor. 

A. Equation for moments 

We are interested in the average quantum probability density {(tP'{x)) = (u^(t)), where 

(M^(r)) — I u^P{u,v,T)dudv . 



From the FPE we obtain a system of equations for the moments 

MkA = {u'v') , (5) 

where k, I = 0,1,2, ... . Substituting u'^v^ in the FPE and integrating over u and v, one obtains the following relation 
for Mkj 

Mk,i = -lijMk+u-i + kMk-u+i + l{l - l)2DMk+2.i-2 + (3lMk+3,i-i , (6) 

where Mk,i with negative indexes are assumed to vanish. We note that only terms with the same parity oi k + I are 
coupled. Since we are interested in M2,o = (w^), we study only the case when this parity is even, namely k + I — 2n 
with 71 = 1,2,.... The sum of the indexes of the moments is 2n, except the last term piM^+s^i-i, where the sum is 
2(71 + 1). This leads to the infinite system of linear equations that can be written in the form 

M = WM , (7) 

where the column vector is M = (M2.0, Mi^i, Afo,2, -^^4.0i M^^i, . . .) and W is the corresponding matrix. The matrix 
elements Wk.i are determined by Eq. ([6]). The solutions of the system of linear equations ([7]) are linear combinations 
of the eigenfunctions at time r 

MA(T)=exp(AT)UA, (8) 

where Ua = Ma {t = 0) is the eigenvector of W corresponding to the eigenvalue A found from the equation 

WVm = XmVm , U„, = Ua,„ . (9) 



The infinite matrix W consists of two parts. The first one is independent of (3 and consists of independent diagonal 
blocks An of size (2n + 1) x (2n + 1). The second one consists of the (3 dependent terms which couple the n-th 
and the {n + l)-th blocks and are located above the (n + 1) block and to the right of the n-th block. Consequently, 
the P dependent terms do not affect the characteristic polynomial, as can be shown by elementary operations on 
determinants. Therefore, the characteristic polynomial of W reduces to a product of the block determinants [20J 

oo 

J|det(A„-A/„)-0, (10) 

n=l 

where /„ is an (2n+ 1) x (2n+ 1) unit matrix. The diagonal block A„ of the infinite matrix W defined in equation ([7]), 
that couples the moments of order 2n with one another, is a band-diagonal square-matrix of size {2n+ 1). The explicit 
form of this matrix is given by: 

/ 2n \ 



A^ri 






2n 







—UJ 





2n-l 





2D 


~2uj 





2n - 2 





6D 


—Slo 


2n-3 



(11) 



(2n-2)(2n-l)i:> -(2n-l)cj 1 

V {2n- l)2nD -2nuj J 

Let us denote by Aniax('T-) the maximal eigenvalue of this matrix. The vector of the moments in this block is M„ = 
(Af2n,0j -^2n-i,i ■ • ■ , M2n-i,i, ■ ■ ■ i -^'^o,2n)- 1^ Appendix A it is proven that for uj = the maximal Lyapunov exponent 
Amax('T-) behaves for large n as 



4 
Then it is argued and verified numerically that also for other values uj ^ it behaves in this way. 



A,„ax(n) ~ pi/3(2,,)4/3 (12) 



III. ASYMPTOTIC GROWTH OF THE MOMENTS 
A. Eigenvalue problem for the moments 

Taking into account Eqs. ^ and ^, we present the solution of Eq. ([7]) as an expansion 

M(r) - Y. Cm(T)U,„ = Y, e^""c„,U„, , (13) 

m ni 

where Cm = Cm{T = 0). Due to the block structure the eigenvectors are characterized by two indexes m = (n^k), 
where n indicates number of block, while k = 1,2,... ,2n-fl counts elements inside each block. Therefore, the 
eigenstates U„i = XJn,k are found from the following algorithm. For the block n — 1 there are 3 eigenvalues Ai_fc 
with corresponding eigenvectors Ui,fc determined by the first block Ai. Therefore, Ui_fc — (Uijt,0), where O is 
an infinite zero vector. For the second block n — 2 there are 5 eigenvalues A2,fc with corresponding eigenvectors 
U2,fc = (Ri,fc, U2,fc, O), where Ri.fc is a 3 dimensional vector, while U2,fc is a 5 dimensional vector, and fc = 1, 2, . . . , 5. 
Here A2.fc and \J2,k are determined from the second diagonal block matrix A2, while Ri,fc is determined by Ai and 
by the corresponding /3-dependent off diagonal block. Continuing this procedure, we obtain 2n + 1 eigenvectors for 
Xn,k in the form 

U„^fe-(R„_i,fe,U„,fc,0), (14) 

where R„_i_fc is a (71^ — 1) dimensional vector determined by n — 1 diagonal and off diagonal blocks of the truncated 
matrix W. 

Summation over m in Eq. (J13p is broken into the sum over the block numbers n e [1, 00) and the sum over indexing 
inside each block I € [0, 2n]. Thus, Eq. (fT5|) reads 

00 2n 

M(^)=m]^".'e^"''"u„,z. (15) 

n=l ;=o 



The vector M consists of the block vectors M„: M = (Mi, M2, . . . , Mji, . . . ), where M„ is a vector of 2n+l elements 
defined in Eq. ^ and corresponds to the moments of the order of 2n. Therefore, the initial vector at r = is 

(Ml, M2, . . . , M„, . . . ) = 5m c«',pU„-,,, . (16) 

n' I' 

A ssume that at the initial point r = 0, the wave function and its derivative are small, of the order of e in units of 
^y\uJ/|3\. Then, at that point the moments scale as M„ ~ e^" with e arbitrary small. As follows from Eq. (|16p. 
one finds Cn,iUn,i ^ Cn^^^ + o(e^") with bounded c„ for any n, as demonstrated in Appendix B. In the linear case 
(/3 = 0) the growth rate of each moment of the order of 2n corresponding to the n-th block is determined by the 
eigenvalue with the largest real part, and we denote it by Aniax('T-) = max(ReA„_;), where the maximum is over the 
2n + 1 eigenvalues corresponding to the nth block, indicated by Z. As shown in Appendix [SI the asymptotic behavior 
of the generalized Lyapunov exponent Amax("-) when n ^ 00 is given by Xmax{n) ~ An'^'^, where A — !^2'^''^D^'^ is a 
constant (see Eq. (fT^ V The leading contribution to the growth of M„(t) in the nonlinear case (/3 7^ 0) is determined 
by the following sum 

M„(r) - ^ c„e2-e'^'"'''^ (17) 

as is clear from Eqs. (fTS jI - fTB]) . 

B. Resummation 

This series in Eq. (I16p has a vanishing radius of convergence and probably has to be interpreted as an asymptotic 
series. It can be used to study the behavior of M„(t) after being resummed. Such a resummation is done with the 
help of the identity 

1 C^°° 
exp(i^2) = -= / dwexp(-u2 + 2Ku) , (18) 

known as the Hubbard-Stratonovich transformation. We can rewrite the above series as follows 

M„(t) - y c^e^™ exp(y^m4/V) = -^ / duexp{~u^)<P{2uVA^) , (19) 

^ V'^ J -00 

m>n 

where the function $ is given by 

*(y)= E^™e'"exp(m2/3y). (20) 

rrt>n 

If the coefficients Cm do not grow too fast (which we shall assume hereafter) the function <&(?/) is well defined at least 
for small values of y. Note that for the resummation of Eq. IT^ it was not crucial that the power of m is 4/3. Such 
a resummation can be performed for any power a < 2 replacing 4/3. 

A more explicit resummation procedure for the r.h.s. of equation p7p can be developed with the help of fractional 
derivatives. First, let us expand the exponential function 

exp(^rm^/3) = ^1 k^ ' ' (21) 

k=0 

Then, writing e^™ = e™'"^ = e""™, we obtain 

Mn(-) = E ^ E cn^m''^'^^"^ . (22) 

A;— m'>n 

We now introduce the Weyl fractional derivative of order g of a function f{z) by the Weyl integral, see e.g. [33 |. 

d'fiz) _ 1 r f{y)dy 



dzi T{-q)J_^{z-yY+^ 



(23) 



where for q > Q the integral should be properly regularized J32|.|33J. and r(— g) is the gamma function. For f{y) ~ e'^'» 
it takes the form 

Substitution of Eq. ^^ in Eq. (gH) with v = m yields 

fc— m>n 

If Cm are bounded, as shown in Appendix B, by some (?„, the sum in Eq. (|25p is bounded by %^^-s-, hence ]VI„ is the 
fractional derivative of some well defined function presented in Eq. (j23p with its regularization. 

Therefore, Eqs. ((25)) and P^ describe the long time asymptotics of the moments. Therefore, Eq. (fT7|) is an 
asymptotic expansion and is a good approximation as long as we sum decreasing terms. The condition is (for 
bounded Cm) 



e^" exp I 



(Am^'^T\ > e^^^+i) exp (A{m + l)^^\) . (26) 

For large m this inequality is In - > ■^ATm^''^ with m > n. Consequently, for the time of the order 

(n)_3 1n(l/e) 

the n-th moment will be dominated by the leading terms and will grow as in the linear case (/3 = 0). 

C. The growth of the second moments 

The second moments are of particular interest for the present work. Their growth for a time that is shorter than 
Tg of Eq. (|27p is dominated by the leading term, namely M2 = C2e^e'^"""=(^''^. Consequently, for r < t^ , 

(u^) = Af2,o = cae'e^— (1)^ . (28) 

This result was verified numerically. Using the analogy between the stationary Schrodinger equation ([1]) and the 
Fokker-Planck equation ([U we identify the generalized Lyapunov exponent Aniax(l) with the growth rate ^ as 

27 = A,nax(l) - lim hm ^^('^'(^)/^') > . (29) 

For T > Tq , non-linear saturation effects become relevant and the full nonlinear theory should be used. 

IV. SUMMARY 

In this work the asymptotic behavior of the generalized Lyapunov exponents of the linear Fokker-Planck equation 
([3]) with [3 — 0, was found analytically and is given by Eq. P^ . The resulting expression for the moments (fT7|) is 

divergent, but it can be resummed in the form Eqs. (|19p or (|25p . Therefore, for short time r < Tq , the first term in 
Eq. p7p provides a good approximation of the moments. In particular, for the second moment this result enabled to 
identify the generalized Lyapunov exponent of the linear system with the asymptotic growth rate for the nonlinear 
one, namely 

27 = Amax(l) ■ 

According to the implementation of the d la Borland method, as outlined in [20|, this is the decay rate of the stationary 
states of the random nonlinear equation ([T]) , showing that it is independent of the nonlinearity /3. 
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APPENDIX A: ASYMPTOTIC BEHAVIOR OF THE GENERALIZED LYAPUNOV EXPONENTS 

In this appendix, Amax("), the maximal eigenvalue of the matrix (jlip will be evaluated. Following [28] we call 
this quantity a generalized Lyapunov exponent. By elementary dimensional analysis [28i |. one finds that the following 
scaling is suggestive, 



^W = D'^'L- (^) • (Al) 



In the long time limit the moments of order 2n grow as e'^""*'=("'*. In [23], Zillmer and Pikovsky studied Xma^in) for 
different values of n: they give exact expressions for n = 2, for n ^ (which corresponds to the usual Lyapunov 
exponent) and consider limiting cases for large values of the dimensionless parameter uj/D^^'^. They also study 
numerically the behavior of Amax(^) as n —> cxd, keeping the dimensionless parameter fixed. They found numerically 
the scaling law: 

Amax(»^) oc n" with a~1.4. (A2) 

The fact that the scaling exponent a is different from 2, implies deviations from Gaussian behavior and breakdown 
of monoscaling; the consequences of this breakdown of single parameter scaling for conductance distribution were 
studied by Schomerus and Titov [29| . We also remark that related studies were also carried out in the context of the 
harmonic oscillator with random frequency [30, IMl ■ 

In this appendix, we study analytically the behavior of the generalized Lyapunov exponents in the limit n ^ oo. 
For the special case of w = 0, we prove the following asymptotic formula: 

A.„ax(n) ~ Ji^l/3(2^)4/3 (A3) 

We shall then argue that this behavior remains valid for any finite value of lo. 

We now outline the proof of the scaling equation (|A3|) for w = by studying the large n behavior of the coefficients 
of the characteristic polynomial P{X) of the matrix An. Recalling that if A„.i, . . . A„.2n+i are the eigenvalues of An 
we have 

2n+l 2n+l 2n+l 

P{X) = n (^ - ^'m) = ^'"+' - ( E ^".0^"" + (E A„,^A„,,)X2"-l + . . . + n Xn,^ ■ (A4) 

i—1 i=l i^j i=l 

We also have X(n) — max/(Re A„^;). The coefficients of the characteristic polynomial P{X) are symmetric functions 
of the eigenvalues Xn,i. Thanks to the Newton formulae, all these coefficients can be written as linear combinations 
of traces of powers of A„. Hence, we have 

P{X) = A:2"+i - Tr(yl„)A:2" + i('[Tr(yl„)]2 - Tr(A2)) AT^"-! + . . . (A5) 

In principle, we can obtain the eigenvalues from the traces of the 2n + 1 powers of A„. In practice, since we are 
interested only in the asymptotic behavior of Ainax('T'), it will be computed from the traces of high powers. From 
equation (jlip we observe that Tr(yl„) = 0. In the case w = 0, we also have Tr(j4^) = 0. More generally, for w = 0, 
we can show that Tr(A^j) ^ only when fc is a multiple of 3. Indeed, writing 

An^d + g (A6) 
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(A7) 



2n{2n - 1) / 



we obtain 



[d + gf 



E 



W^fcl,fc2 



(A8) 



where Wk-i,k2 is a product of k factors with fci factors equal to d and fc2 factors equal to g, with fci + fc2 = fc (there are 
2'^ such terms because the matrices d and g do not commute). We remark that d is an upper-diagonal band matrix and 
its non vanishing terms are all on the band located at level +1 above the diagonal. Similarly g is a lower-triangular 
band matrix and its non vanishing terms are all on the band located 2 levels below the diagonal. Therefore a product 
of fci matrices d and fc2 matrices g will have a non-zero diagonal term only if fci = 2fc2 i.e. if fc = 3fc2; hence, fc is a 
multiple of 3. For example, we have 



2n-l 



TT{Ai) ^TT{ddg + dgd + gdd)^3Ti{d-'g) = 3D > ;(/-|-l)(2n-/-|-l)(2n-0 =i 3i:'(2n)^ / x\l-xydx^ D 



1=1 



16n^ 



(A9) 

More generally, the terms that contribute to Ti-{Aj^) are obtained by taking the product of 2fc factors d and k factors 
g written in all possible orders (there are AkV.ki such terms); 



Tr(ylf ) = Y^ TT{W2k,k) ^ d^'^g'' + 'permuted terms' 



(AlO) 



In particular, we have 

Tr(d2'=gfe) 



k sr^ 



D 



D'=(2n 



l{l + l){l + 2)..{l + 2k 



Ak+l 



„2fe 



(1 



\2k 



l){2n -l + l)i2n- l)...{2n ~ I ~ 2k + 2) 
,4,+i(2fc)!(2fc)! 



dx = D''{2n) 



{Ak + iy. 



(All) 



The trace of any term W2k,k is given by the same expression at the leading order: indeed, the elements of the matrix 
dg are of order tt' whereas those of the commutator [d,g] are of order n^. The reason is that both dg and gd are 
triangular with all nonvanishing matrix elements one level below the diagonal, of the form (dg)i^i-i and {gd)i^i-i. In 
the center of the matrix, Z « n, a generic term is of the order of n^ + O(n^). Therefore both dg and gd are dominated 
by n^, while [d,g], given by the difference of two such terms is dominated by the 0{n^) corrections. We know that any 
W2k,k differs from d^^g'' by a finite number of commutators. Therefore, Tr(W2fe,fc) = Tr{d^''g'^)+ subleading terms. 
We thus have 



Tr(Af) 



D''{2n 



i4fe+l 



(3fc)!(2A:)! 



(4fc-fl) (4/c)!fc! ■ 
Rewriting this trace in terms of the eigenvalues of the matrix An, we obtain 

2«+l ^fc/n„^4fc+l (3fc)!(2/c)! 



Besides, from equation (jA12|) . we deduce that 



3fc 

'n,2 



D''{2n) 
(4fc-f 1) {Aky.kl 



(A12) 



(A13) 



Tr(A^) 

~t^{aW) 



k~l 



(A14) 



2n 


Amax(n) 


4A„ax(n) 


10 


15 


1.077 


40 


101 


1.015 


80 


257 


1.006 


100 


347 


1.003 


120 


443 


1.0021 


140 


544 


1.0022 


160 


650 


1.0018 


180 


761 


1.0016 


200 


876 


1.0013 


300 


1504.7 


1.0010 



TABLE I: Behavior of the dominant eigenvalue in the case uj = Q, D = 1. 



2n 


Amax(".) 


3D 
4A 


2n)''/3 


10 


13.2 




1.223 


40 


95.9 




1.072 


80 


248 




1.043 


100 


336 




1.035 


120 


431 




1.031 


140 


530 




1.028 


160 


635 




1.026 


180 


745 




1.023 



TABLE II: Behavior of the dominant eigenvalue in the case lj = 1, D — 1. 



This equation shows that the trace is not dominated by the largest eigenvalue alone: otherwise this 'participation' 
ratio would be of order one. Rather, a finite fraction p of the eigenvalues has a scaling behavior similar to that of 
Amax('^) and we can write: Tr(^^'^) ~ pnXma.xin)^'' ■ Thus, we deduce from equation (|A13p the following behavior 



3, Jjn2n)^^+i (3fc)!(2fc)! 
Finally, using Stirling formula, we obtain for large fc, assuming pn does not vary strongly with n 



A„,„,.Z,V3p„,.3((3|ffl)./»^3^,,,^^,, 



(A15) 



(A16) 



This ends the proof of equation (|A3|1 or (|12p in the case uj = 0. In Table I the highest eigenvalue of the matrix yl„ is 
computed numerically for various values of 2n up to 2n = 300 (for a; = and D = I). The asymptotic scaling given 
in equation (|A3p is well satisfied. 

When Lu is different from 0, we can still write A„ as a sum of two matrices as in equation (|A6P : the upper-diagonal 
matrix d remains the same as in equation ()A7p but the lower-triangular part g' is now given by the previous g plus 
a band-diagonal matrix (72 (w) (with a band-diagonal at level -1) which contains terms proportional to uj. However, 
in the large n limit (and keeping the value of w fixed) the matrix elements of g are much bigger than those of g2{^)'- 
hence (7' ~ g up to subdominant contributions and the large n scaling of the maximal eigenvalue is insensitive to lo 
at leading order, therefore equation (|A3p . derived for w = 0, remains true for finite values of lo. In Tables II and III 
we give numerical results uj = 1 and uj = —1 (taking D — I for both cases). The scaling behavior, proportional to 
(2n)''/^, is well satisfied and it seems also that the prefactor 3/4 remains correct. 



2n 


Amax(n) 


3D 

4A 




10 


18 




0.900 


40 


107.5 




0.952 


80 


267 




0.968 


100 


358 




0.972 


120 


455 




0.975 


140 


557.5 




0.977 


160 


666 




0.980 


180 


777 




0.981 



TABLE III: Behavior of the dominant eigenvalue in the case u) — —1, D = 1. 



APPENDIX B: SOUNDNESS OF THE c„ 

Let us assume an expansion (of the initial moments (116p l 

oo 



(Bl) 



m=0 



for the l-th moment of order n as defined after the matrix pip . The 6's are bound (by construction). Let us expand 
the Cn.i of Eq. (fT6|) in powers of e as 



CnA 



EC 



(m') 2m' 



(B2) 



Since only even moments are considered, the powers of e are even. Now we write Eq. p6p in the form 

M„.( = ^c„^i.C/^?;9, (B3) 



n'V 



(n.l) 



where we equate the components of the vectors, now we use the property (|14p of the eigenvectors C^„?';/ , namely 
Uz?], ~ for n > n' , and solve the equation order by order in e. 



Take first n ~ 1, and leading order in e 



6^'^- Vc^'^C/^'''^ /-012 
1,' ~ Z^ i,i' 1,'' t — U, 1,Z. 



(B4) 



These are three linear equations for the ci^p, since the eigenvectors U^ ;', are given. 
Next take n = 2, 



^S = E€^<^ + E€^<^ 



l'=0 
2 



4 



bS-T.^'US'' + T.^'US' 



(B5) 
(B6) 



l'=0 



l'=0 



_(2) _(2) 

These are eight equations for the C2 // , €[ f, . 



This process can be continued to any order in e. The independence of the eigenvectors implies that finite solutions 
] can be obtained by Kramer's rule. Note, that b, ' = for n > m' therefore, also c , ^ = for 71 > m'. 



Therefore, 



Cn,Z ^^ Cn^ 



(B7) 



for small e. 



This shows that the c„ are bounded for any finite n, but it does not imply the existence of a uniform bound for all 
n and m. 
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